%% 
% 'a' needs to be Nx2 array, where first column is temperature, second is
% growthrate
load('gr_data_all_20230413_v2.mat');
close all
a = gr_data_all;
T = unique(a(:,1));

MGcolor = [235, 147, 52]/255;

G = [];
Gerr = [];
for k=1:6
    g = reshape(a(:,k+1),[length(a)/length(T) length(T)]);
    gmean = nanmean(g,1)';
    gstd = nanstd(g,1)';

    G = [G gmean];
    Gerr = [Gerr gstd];

end


%% Plot data
T1plot = 16;
T2plot = 50;

clrs = distinguishable_colors(6);
clrs = cbrewer('qual', 'Paired',6);
strainnames = {'MG1655', 'BW25113', 'CS109', 'BL21', 'dnaK', 'uup'};

ms = 6;
figure;

toplot = [1];

for k=toplot
    g_plot = G(:,k);
    gerr_plot = Gerr(:,k);
    errorbar(T,g_plot,gerr_plot,'ko', 'markersize', ms, 'markerfacecolor', MGcolor);
    hold on;
    
end

set(gcf, "Position", [0 0 400 300]);
set(gca, 'FontSize', 20);
xlabel('Temperature (°C)');
ylabel('Growth rate (1/h)');
box off;
legend(strainnames{toplot})
xlim([T1plot T2plot])

%% Plot Arrhenius data

Tmin = 25;
Tmax = 37;
clc;
Ea_tot = [];
Ea_err_tot = [];


figure;
for k=toplot
    g_plot = G(:,k);
    gerr_plot = Gerr(:,k);
    errorbar(1./(T+273),log(g_plot), gerr_plot./g_plot, 'ko', 'markersize', ms, 'markerfacecolor', MGcolor);
    hold on;

    Tfit = T(T<=Tmax & T>=Tmin);
    gfit = g_plot(T<=Tmax & T>=Tmin);
    gfit_err = gerr_plot(T<=Tmax & T>=Tmin);
    
    Tfit = Tfit(~isnan(gfit));
    gfit_err = gfit_err(~isnan(gfit));
    gfit = gfit(~isnan(gfit));
    

    % Linear fit
    [f1,f2] = fit(1./(Tfit + 273.15), log(gfit), 'poly1');

    % Weighted linear fit
    ft = fittype('poly1');
    cf = fit(1./(Tfit + 273.15), log(gfit),ft,'Weight',gfit_err./gfit);
    
    % Weighted linear fit using fitnlm
    Arrhenius_fit = @(b,x) b(1) + b(2).*x;
    x0 = [30 -8000];
    f = fitnlm(1./(Tfit + 273.15), log(gfit),Arrhenius_fit,x0, 'Weights', gfit_err./gfit);

    Tplot = linspace(25,37,20);
    hold on;
    plot(1./(Tplot + 273.15), f.feval(1./(Tplot + 273.15)),'k', 'linewidth', 1)

    Ea_tot = [Ea_tot; -f.Coefficients.Estimate(2)/298/1.688];
    Ea_err_tot = [Ea_err_tot; f.Coefficients.SE(2)/298/1.688];
end
set(gcf, "Position", [0 0 400 300]);
set(gca, 'FontSize', 20);
xlabel('1000/T (1/K)');
ylabel('Growth rate (1/h)');
xlim([1/(T2plot + 273.15) 1/(T1plot + 273.15)])
box off;